/*  sdfatan2.cpp  by K.Tsuru  */
// function ID = 3110  DARDIX
/*******************************************************
A bug for x = 0 was removed in version 2.14.
SDouble class
inverse trigonometric function arctan(y/x)
SDouble version of atan2(double y,double x)
Return an azimuth angle of the point (x,y) on the xy-plane
within the region (-pi,pi].
********************************************************/
#ifndef SN_H
#include "sn.h"
#endif

#if 1
SDouble Atan2(const SDouble& y, const SDouble& x){
    SDouble r;
    if(x.Sign() == 0){
        r = MPi2(); // r = pi/2
        if(y.Sign() < 0) r.ChangeSign(); // r = -pi/2
        if(y.Sign() == 0) r.SetError(x.DOMAIN_ERR, "Atan2", -3110);
        return r;
    }
//x != 0
    if(y.Sign() == 0){
        if(x.Sign() > 0) return 0.0;
        return Pi();
    }
    r = y/x;
    r = Atan(r);
    if(x.Sign() > 0);  // r = Atan(r);  1st or 4th
    else if(r.Sign() < 0) r += Pi();    // x<0, y >0, r = Pi() + Atan(r);   2nd
    else r -= Pi();     // r = Atan(r) - Pi(); 3rd or y = r = 0
    return r;
}
#else
static SDouble pi= Pi();
SDouble Atan2(const SDouble& y, const SDouble& x){
    SDouble r;
    if(x.Sign() == 0){
        r = MPi2(); // r = pi/2
        if(y.Sign() < 0) r.ChangeSign(); // r = -pi/2
        if(y.Sign() == 0) r.SetError(x.DOMAIN_ERR, "Atan2", -3110);
        return r;
    }
//x != 0
    if(y.Sign() == 0){
        if(x.Sign() > 0) return 0.0;
        return Pi();
    }
    r = y/x;
    r = Atan(r);
    if(x.Sign() > 0);  // r = Atan(r);  1st or 4th
    else if(r.Sign() < 0) r += pi;  // x<0, y >0, r = Pi() + Atan(r);   2nd
    else r -= pi;       // r = Atan(r) - Pi(); 3rd or y = r = 0
    return r;
}
#endif